Initialize printing:


In [1]:
%pylab inline
from sympy.interactive.printing import init_printing
init_printing()


Welcome to pylab, a matplotlib-based Python environment [backend: module://IPython.zmq.pylab.backend_inline].
For more information, type 'help(pylab)'.

Poisson Equation

Radial Poisson equation is $$ \phi''(r)+{2\over r} \phi'(r) = -4\pi\rho(r). $$ Alternatively, this can also be written as: $$ {1\over r}(r \phi(r))'' = -4\pi\rho(r). $$

Example I

Positive Gaussian charge density:


In [2]:
from sympy import var, exp, pi, sqrt, integrate, refine, Q, oo, Symbol, DiracDelta
var("r alpha Z")
alpha = Symbol("alpha", positive=True)
rho = Z * (alpha / sqrt(pi))**3 * exp(-alpha**2 * r**2)
rho


Out[2]:
$$\frac{Z \alpha^{3} e^{- \alpha^{2} r^{2}}}{\pi^{\frac{3}{2}}}$$

The total charge is $Z$:


In [3]:
integrate(4*pi*rho*r**2, (r, 0, oo))


Out[3]:
$$Z$$

Solve for $\phi(r)$ from the Poisson equation:


In [4]:
phi = integrate(-4*pi*rho*r, r, r)/r
phi


Out[4]:
$$\frac{Z \operatorname{erf}{\left (\alpha r \right )}}{r}$$

Tell SymPy that $\alpha$ is positive:


In [5]:
phi = refine(phi, Q.positive(alpha))
phi


Out[5]:
$$\frac{Z \operatorname{erf}{\left (\alpha r \right )}}{r}$$

Plot the charge:


In [6]:
from sympy import plot

In [7]:
plot(rho.subs({Z: 1, alpha: 0.1}), (r, 0, 100));


Plot the potential


In [8]:
plot(phi.subs({Z: 1, alpha: 0.1}), 1/r, (r, 0, 100), ylim=[0, 0.12]);



In [8]: